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Abstract 

The Galerkin method is used to derive a reahstic model of plane 
Couette flow in terms of partial differential equations governing the 
space-time dependence of the amplitude of a few cross-stream modes. 
Numerical simulations show that it reproduces the globally sub-critical 
behavior typical of this flow. In particular, the statistics of turbulent 
transients at decay from turbulent to laminar flow displays striking 
similarities with experimental findings. 

PACS: 47.27.Cn, 47.60.+i 

1 Introduction 

The transition to turbulence in shear flows close to a solid wall is far from 
being completely understood. This situation arises because linear stability 
analysis is less fruitful for systems in which transient algebraic energy growth 
may become relevant than when exponentially growing modes are present pQ . 
It seems indeed easier to understand the case of super-critical instabilities 
(especially in closed-flow configurations such as Rayleigh-Benard convection 
and the Taylor-Couette instability for which the classical tools of weakly 
nonlinear analysis are available) than the case of discontinuous transitions 
marked by a competition between solutions derived from linear theory and 
other fully nonlinear solutions. An interesting review of the main issues 
related to this problem can be found in [2]. A classical example of globally 
sub- critical transition is plane Poiseuille flow — the flow between two fixed 
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parallel plates driven by a pressure gradient — which is linearly unstable only 
beyond some high Reynolds number [3] while turbulent spots, i.e. patches 
of turbulent flow scattered amidst laminar flow and separated from it by well 
deflned fronts, can develop for R as low as about -Rc/4 [1]. The situation is 
even worse for plane Couette flow (pCf in the following), the simple shear 
flow driven by two plates moving parallel to one another, which is stable for 
all R, i.e. Rc = oo [S], but experiences a direct transition to turbulence at 
moderate values of i? [HI [71 [El [9] . Linear theory, modal and non-modal [l] , 
arguably indicates ways to escape from the laminar regime but is of limited 
help to understand the main problem alluded to above, namely the nature of 
the nontrivial nonlinear state and the coexistence in physical space of these 
different solutions. 

Direct numerical simulations of the Navier-Stokes equations (DNS) have 
been intensively used to identify elementary processes involved in the suste- 
nance of the nontrivial turbulent state, see e.g. An interesting approach 
was the reduction to a minimal flow unit (MFU) introduced by Jimenez and 
Moin [TT], in which the dimensions of the simulation domain were decreased 
down to sizes for which turbulence is just maintained. The resulting dynam- 
ics can then be understood as a chaotic evolution of large coherent structures 
driving smaller eddies in a more or less stochastic way (though strict sepa- 
ration of scales is not legitimate). From observations of near- wall turbulent 
flow [TUj, Waleffe was able to derive a differential model involving the ampli- 
tude of a few modes associated to such coherent structures [12] , in the spirit 
of a Galerkin approach that was made more systematic by Eckhardt and 
co-workers [HI [15] using Fourier modes or by Holmes and co-workers [16] 
using empirical modes extracted from DNS by proper orthogonal decompo- 
sition. 

Such models indeed help one to illustrate some of the mechanisms in- 
volved in the sustenance of turbulence or the competition between laminar 
and turbulent states in phase space. However, they are intrinsically and ex- 
plicitly low-dimensional, which situates the problem within the context of 
discrete dynamical systems appropriate to temporal chaos, e.g. in discussing 
transients in terms of fractal border of the laminar-flow basin of attraction 
in phase space [T7] . Therefore, they do not allow one to approach such prob- 
lems as spot propagation that require full reference to the physical space. 
When similar questions were posed for convection, i.e. pattern formation 
and spatio-temporal chaos, it was found particularly valuable to pass from 
ordinary differential systems that can only deal with temporal behavior, such 
as the Lorenz model, to partial differential equations that also account for 
spatial dependence, e.g. the Swift-Hohenberg model pjj and its numerous 
variants. The present paper is devoted to the derivation and use of such a 
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closed set of partial differential equations appropriate to the pCf. This case is 
particularly interesting since, in contrast with other wall-bounded flows such 
as Poiseuille flow or the laminar boundary layer (Blasius) flow, it completely 
lacks linear instability and the absence of downstream advection makes it 
practical to observe the long term dynamics of spots, a crucial element of the 
sub-critical transition to turbulence. 

Contrasting with the regime of developed turbulence taking place at large 
R, which requires a refined space-time resolution able to account for a cas- 
cade toward small scales, the transitional regime around Rg involves struc- 
tures that appear to be coherent and large, i.e. occupy the full gap between 
the plates [19]. In turn, a model in terms of a few well-chosen modes with 
low wall-normal resolution should be sufficient for describing spots and un- 
derstanding their dynamics. The amplitudes of these modes would then be 
taken as functions of the in-plane coordinates. This is the reason why we 
speak of 2. 5- dimensional models. In this terminology, the number 2 stands 
for the full in-plane space dependence (x, z) and the suffix .5 suggestively 
expresses that the dependence on the third, cross-stream, coordinate y is 
only partly taken into account through low order truncation. This strategy 
has been developed previously for easy-to-treat but unphysical stress-free 
boundary conditions at the plates bounding the flow [20] yielding interesting 
results but at unrealistically low Reynolds numbers [21], [22]. Here we focus 
on realistic no-slip boundary conditions. In the next section, we recall the 
primitive equations and give a hint on how the Galerkin method is devel- 
oped. Next we turn to the explicit derivation of the model (the no-slip and 
stress- free models are compared at a formal level in §2.3p . The numerical im- 
plementation is described in ^ Our main results are presented in Section HI 
essentially bearing on the globally sub-critical character of the transition to 
turbulence ( §4.ip . The behavior of associated transients is then studied by 
experiments where the system is quenched from a turbulent state at large R 
towards ever decreasing values of R ( §4.4p . These results are discussed and 
some conclusions are drawn in Section [5l 

2 The 2.5D models of 3D transitional flow 
2.1 Primitive perturbation equations 

The Navier-Stokes equation and continuity condition for an incompressible 
flow read: 



dtv -\- V ■ Vv 
V- V 



-Vp + z/VV + f , 

0, 




(2) 
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where v = {u,v,w), p is the pressure and u is the kinematic viscosity. Fur- 
ther, denotes the three-dimensional Laplacian and the term f accounts 
for a driving bulk force, if necessary. 

Two types of boundary conditions can be considered, either the easy-to- 
handle but unrealistic stress-free (sf) conditions, or the natural and realistic 
no-slip (ns) conditions. In the sf case, the fluid is supposed to slip on the walls 
so that the flow has to be maintained by a fictitious bulk force f oc sm{(3y/h) 
with P = 7r/2, further adjusted to produce the basic velocity profile U^^{y) = 
Up sm[f3y / h) . In contrast, with ns boundary conditions, the motion of the 
plates is able to drive the fiow, which (when laminar) results in the linear 
velocity profile U{^^{y) = U^y/h, where 2f/p is the relative speed of the plates, 
and 2h is the gap between them. 

In the following we use dimensionless quantities. According to common 
usage, lengths and speeds are scaled with /i, and ?7p, respectively. In the 
ns case, the basic velocity profile thus reads ^7^'' = y for y E [—1,1], and 
the Reynolds number is R = Uph/u. The models are further written for 
the perturbation {u' , v' , w' , p') to the laminar basic fiow that defines the x 
direction, Vb = U^Si, with f/b = or U{^^, i.e. u = U^^y) + u', v = v', 
w = w', p = p', so that ( fT|2i) then read: 

dtu + u'dxU + v'dyU + w'dzu' 

= -dj - U^d^u' - v'j-U^ + R-^v\' , (3) 
dtv' + u'dxv' + v'dyv' + w'dzv' 

= -dyp'-U^d,v' + R-'vW (4) 
dtw' + u'dxw' + v'dyw' + w'dzw' 

= -dzp' - U^dxw' + R-^V'^w' , (5) 
= dxu' + dyv' + dzw' . (6) 

The Galerkin method is a special case of a weighted residual method [23] . 
It consists here in forcing the separation of in-plane and wall-normal coordi- 
nates by expanding the perturbations {u' , v' , w' , p') onto a complete basis of 
?/-dependent orthogonal functions satisfying the BCs with amplitudes depen- 
dent on {x,z,t). The equations of motion are then projected onto the same 
functional basis, using the canonical scalar product. The main modeling step 
is then performed when truncating these expansions at a low order and keep- 
ing the corresponding number of residuals in order to get a consistent and 
closed system governing the amphtudes retained. 
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Figure 1: Profiles of tlie basis functions used for the perturbation to the 
base flow in the derivation of the ns model (solid lines), compared to their 
counterparts for the sf model (dashed-dotted hues) using the same mean- 
square normalization: drift flow B{1 — y'^) vs. l/-\/2, first odd component 
Cy{l — y^), vs. sin(7r?//2), and first even component A^l — y"^)"^ vs. cos(7ry/2). 

2.2 The no-slip model at lowest order 

The free-slip model introduced in [20] and used more extensively in [22j was 
developed by taking advantage of the fact that trigonometric functions form a 
closed set under differentiation and multiplication. Perturbations were taken 
of the form: v' = Vi cos{(3y) ancfl {u' , w'} = ^/i/2{Uo, Wo}+{Ui, Wi} sm{(3y) 
with /3 = 7t/2 and the derivation of the model was straightforward. Numerical 
simulations have shown that it displays the expected globally sub-critical 
transition to turbulence but at unrealistically low Reynolds numbers [21] . 
This can be attributed to the underestimation of viscous dissipation effects 
implied by the choice of sf boundary conditions, in the same way as for 
convection where the sf threshold is about 2/5 that for ns conditions (277r^/4 
vs. 1708). The derivation of a model for the realistic case by a similar 
Galerkin approach is harder because simple trigonometric expansions do not 
satisfy the simultaneous conditions: 

v'{y = ±1) = dyv'iy = ±1) = (7) 

obtained by combining the continuity equation (El) with the conditions 

u'{y = ±l)=w'{y = ±l) = 0. (8) 

To overcome this difficulty we take a basis consisting of polynomials in the 
wall-normal coordinate y, which also forms a family of functions closed under 

^Here the amplitude of the j/-independent function is normalized differently than in [201 
I22j in order to ease forthcoming comparisons. 
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multiplication and differentiation. A strict Galerkin approach is followed 
Chap. 11, §10], with projections defined using the canonical scalar product: 

{f,9) = £'f{y)9{y) dy. 

From conditions ([7]) it is clear that the expansion for v' must be taken of 
the form [23 P- HO]: 

v' = ii-yyj2Vn+My) (9) 

n 

where the Rn{y), n = 0,1,... are polynomials of increasing degree n. The 
basis for u' and w' is then obtained by requiring that ([6]) be verified for 
each n and that the basis be orthonormal. It can be checked that these two 
conditions can be satisfied simultaneously and that u' and w' expands as 

W'} = (1 - y') ^ {Un, Wr,} SrM (10) 
n 

where the expression of Sniy) derives from that of Rn{y) using the rules stated 
above. The Galerkin procedure isolates the dependence of the velocity field 
on the cross-stream coordinate y and, in expansions (MTOl) . Um, Kn, Wm are 
just functions of the in-plane coordinates x, z, and time t. 

It is quite simple to restrict consistently to the lowest possible order, 
i.e. keeping only functions associated to Uq, Wq, Ui, Vi, and Wi as for 
the lowest order stress-free model. This is due to the parity properties of 
the functions involved that automatically guarantee the orthogonality of the 
different contributions to the velocity field. Our expansion will thus be based 
on 

v' = V,Ail-yr, (11) 
as the equivalent of Vi cos{j3y) and accordingly: 

K, w'} = {Uo, Wo}B{l - y^) + {f/i, W^}Cy{l - y^) , (12) 

in replacement of -^{Uq, Wq} + {t/i, Wi} sm{[3y). We need not worry about 
expanding the pressure since, at each stage of the Galerkin projection, pres- 
sure contributions — here Pq and Pi — introduce themselves as Lagrange 
multipliers serving to fulfill the successive projections of the continuity equa- 
tion ([6]), here just the projections on the lowest order even and odd velocity 
basis functions fll2l) . Physically, the contribution to the flow identified by 
(?7o, Wq) corresponds to the streak component induced by the lift-up mech- 
anism acting through the component Vi of the contribution (f/i, Vi, W^i), a 



6 



part of which corresponds to streamwise vortices, flow structures that play 
an important role in the self-sustainent of wall-turbulence. Going beyond 
lowest consistent truncation order is tractable but somewhat tedious and 
we believe that the lowest order model contains all the large scale features 
present in the experiment at a qualitative level and can account for them at 
a semi-quantitative level. 

Throughout the derivation, integrals to be computed are of the form: 



2k + n + l 



where the are the binomial coefficients. For example, the normalization 
constants are given by: = l/2Jo,4 = 315/256, B"^ = l/2Jo,2 = 15/16, and 
C2 = 1/2J2,2 = 105/16. 

Amplitudes introduced in flllfl2p are all functions of x, z, and t. The 
drift-flow component {Uq, Wq) now has a plane Poiseuille profile, in close 
correspondence with what happens in the Rayleigh-Benard first no- 

ticed by Siggia and Zippelius [26] and subsequently exploited to derive the 
generalized Swift-Hohenberg model described in [27]. Figure [T] compares the 
profiles of the basis functions for the free-slip and no-slip cases used at lowest 
consistency order. 

According to the prescriptions of the Galerkin method, we insert the 
assumed expansions ( fTT1[T2|) in the continuity equation ([6]), then we multiply 
it by B{1 — y^) and integrate over the gap, which extracts its even part: 

d,U, + d,W^ = Q. (13) 

In the same way, multiplying ^ by Cy{l — y"^) and integrating extracts the 
odd part that reads: 

d,Ui + d-Mi = f3Vi, (14) 

where the quantity (3 = -\/3 = has the same order of magnitude as the 
one appearing in the stress-free model = n/2 (see §2.3p . 

Doing appropriate manipulations on ([3]) and ([5]), we obtain the even parts 

as 

dtUo + = -d^Po - aAUi - aaVi 

+ i?-^(A-7o)f/o, (15) 
dtWo + Nwo = -d,Po-a,d^,W, 

+ R~\A-^o)Wo, (16) 
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with ai = 1/v^, aa = ^27/28, 70 = 5/2, and 

Nu, = a,{Uod,Uo + WodMo) 

+ a2{Uid^Ui + Wid^Ui + /5Vi[/i) , (17) 

Nw, = aiiUodM^ + WodMo) 

+ a^iUAW^ + W^ia.m + P'ViW{) , (18) 

with ai = 3^15/14, ^2 = ^15/6, (3' = 3^3/2 = §/?. In ( JTBllTHll and below 
A denotes the two-dimensional Laplacian c}^,.^,. + dzz- 
In the same way, the odd part of (I3l) and ([5]) read: 

5tt/i + iVc;, = -d,P,-aAUo 

+ i?-i(A-7i)f/i, (19) 
a^H^i + iV^y, = -dzP^-a^d^Wo 

+ R~' (A - 71) W, , (20) 

with 7i = 21/2 and 

Nu, = a2{UodMi + Uid,Uo 

+ WodzU, + W^Uo - /5'Vit/o) , (21) 

+ WodzWi + WidzWo - (3"ViWo) , (22) 

with (3" = V3/2 = 

Finally, the equation for Vi is obtained in the same way by projecting (jlj) 
onto (fTTl) . which yields: 

W + iVy, = -/5Pi + i?-^(A-70V^i, (23) 
Nv, = asiUod,V, + WodzV,), (24) 

with 7; = /?2 and 03 = 5^15/22. 

Notice that coefficients aj behave as constants in the model, while the (5s 
are analogous to wave- vectors accounting for cross-stream differentiation ap- 
plied to the corresponding terms. In the same spirit one could remark that, 
in f|T5|) . a2 = P'ai = |/3ai. Another point worth mentioning is that, by con- 
struction, the nonlinearities in the model preserves the energy conservation 
properties of the advection term v- Vv. This fact is indeed important because 
the nonlinear evolution is formally governed by a quadratic expression which 
may produce singularities in a finite time. Energy conservation can be verified 
by computing the evolution of = Eq + Ei = liU^ + W^) + l{Uf + + W^) . 
This tedious task can be by-passed by taking advantage of the identity: 
V • Vv = V (|v^^ + (jj ■ V, with (jj = (Vv)* — Vv, where (Vv)* is the 
transpose of tensor Vv, which leads to the rewriting of the advection terms 
given in the appendix. 
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Table 1: Comparison of coefficients for the sf and ns models at lowest order. 
Definitions are those introduced for the no-slip model and further identified 
term-by-term in the stress-free model (the factors a/2 all stem from the nor- 
malization convention for the y-independent mode associated to amplitudes 
{Uq, Wq}, see Figure [Hand the text relative to the sf basis). 



coefficient 




f3' 


(3" 




a2 


stress-free (sf) 


tt/2 







1/V2 




no-slip (ns) 


Vs 


1/5 




1/V7 


3 „ ns /Q ns 


sf/ns 


0.907 


0.605 





1.871 


1.131 



coefficient 


ai 


a2 


as 


7o 


71 




stress-free (sf) 




1/V2 


1/V2 









no-slip (ns) 


3^/15/14 


VTE/6 


5VTE/22 


5/2 


21/2 




sf/ns 


0.852 


1.095 


0.803 





0.235 


0.822 



2.3 Comparison of stress- free and no-slip models 

The no-slip model derived above appears to be slightly more general than 
the stress- free model initially presented in [20] . The two models have exactly 
the same structure, apart from the term multiplied by jS" which is absent in 
the stress-free case. Furthermore, when normalized in the same way, their 
coefficients are comparable, while the main difference lies, as expected, in the 
terms associated to viscous dissipation. Let us examine Table [T] in detail. 

As already said, values found for the quantity P in ( fT^ . which measures 
the typical order of magnitude of gradients in the wall-normal direction, are 
close to each other since /3'^^ = ^3 ^ 1.732 and /3'^ = tt/2 ^ 1.571. 

Coefficients ai^2 and ai_5 appearing in the other equations are also not 
very different in the two models. From the sf/ns entries in the table, one can 
see that ratios are all within a factor less than two, with the exception of 
f5" which is zero in the stress-free case owing to the special relations existing 
between trigonometric lines. As far as non-normal effects and nonlinear in- 
teractions are concerned, one may think that the two models are close to each 
other and that the dynamics that they generate will be robust (the effects 
of the differences could be studied in detail by considering fictitious alter- 
nate models with arbitrary coefficients, i.e. not derived from any systematic 
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Galerkin approach). 

The most important difference between the two models shows up in the 
expression of the viscous terms. Whereas in the stress-free model, the drift 
velocity component {Uq, Wq} can relax only through in-plane modulations via 
the terms involving A, an additional damping is observed in the no-slip case 
with coefficient 5/2 in fll5fl6l) independent of the in-plane dependence of that 
flow component. In the same way the damping of components {f/i,iyi} is 
much stronger in the no-slip case with coefficient 71 = 21/2 in ( fT9|20l) . than 
in the stress-free case with coefficient (vr/2)^ ^ 2.41. The ns wall-normal 
viscous timescale is thus four times shorter than the sf one. Accordingly, 
if turbulence can be sustained in the stress-free model at some value R, 
already independently of other sources of damping and while the mechanisms 
are expected to retain most of their intensity as discussed in the previous 
paragraph, a similar situation is expected for about 4i? in the no-slip model. 
Of course the argument is very rough but it shows that one goes in the right 
direction since the stress-free model is known to underestimate thresholds by 
a large factor ^T\ . 

Another important feature of the no-slip model is that an (x, z)- 
independent component of Ui can be created as the flow evolves, so that the 
model already contains a mean flow correction to the base profile Ub{y) = y 
even when truncated at lowest order. In contrast, a similar correction for- 
mally appears only at a much higher order in the stress-free model, as pointed 
out in [15J. In the latter case, velocity profiles attached to Ui and Wi are 
indeed the same as that of the base flow and accordingly do not modify its 
shape though an equivalent mean contribution also exists. This feature is 
further examined in the next subsection. 

2.4 Mean flow correction 

Let us consider the average value 



of the streamwise perturbation velocity component Ui over the domain V 
at the boundary of which perturbations cancel, or else periodic boundary 
conditions apply. The equation governing Ui(t) is easily obtained by av- 
eraging (fT9|) over the domain. On its right hand side, the only term that 
remains in the average is the one that is not the gradient of some quan- 
tity, i.e. —■jiR'^Ui, accounting for viscous dissipation. The evaluation of 
Iv Nu^dxdz leads to Nu-^ = a2{l3 + I3")V~^ /j, UoVidxdz. This quantity — the 
space- independent contribution to the streamwise component of the Reynolds 



V 
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stress — is easily obtained after some simple algebra using the continuity con- 
ditions fll3|14l) . Gathering the results we get the evolution equation for the 
mean flow correction as: 




\ = a2{P + P")UoV,-^,R-'U^. 



(25) 



As will be illustrated later (cf. Figure HI), a fully turbulent regime at steady 
state will display a time-independent non-zero — in fact negative — correc- 
tion, arising from the compensation of the two terms on the right hand side. 

The spanwise component of the correction Wi is governed by the same 
equation as ( l25l) provided that we make the general change U ^ W. However 
Wi is expected to cancel on average for a turbulent regime at steady state as 
a consequence of the z —z symmetry of the original problem. The full set 
of perturbation equations is indeed symmetric in the interchanges x ^ z and 
U W except for the presence of the term —02^1 accounting for lift-up in 
(|T5l) . which has no equivalent in its spanwise analogue ( JT6l) . This is the only 
term that singles out the streamwise direction, subsequently introducing the 
difference between Ui and Wi. A non-zero Wi could then only result from 
an instability breaking that symmetryjl 

Whereas it is clear that no such correction exists for Vi since Vi obviously 
cancels all the time (no net flux across the layer for impermeable walls), we 
must also consider averaging the fields Uq, Wq. The equation for Uq reads: 



A similar equation can be obtained for Wq through the change U 1— > W. The 
source terms for Uq and Wq are however expected to be small since, from 
the continuity condition IHM . Ui and Wi are spatially out of phase with Vi. 
These terms can even be expected to cancel statistically since otherwise they 
would generate net contributions to the transfer of momentum in the x and z 
directions and correlative finite dissipation terms —^R'^Uq and —^R'^Wq. 
As long as the global symmetries of the problem are not broken, net fluxes 
U and Wq are forbidden by the general x 1— > —x and z t— s> —z symmetries in 
the same way as Wi 7^ is forbidden by the z —z symmetry (in contrast 
with Ui which has its own source term). 

It can be noticed that, in the stress-free case, mode '0' has no source term 
for accidental reasons linked to trigonometric relations between base func- 
tions (i.e. (3 = 13'). Furthermore, that mode has no cross-stream dependence 
so that it cannot be ironed out by viscous dissipation. Accordingly, Uq and 

■^Such a symmetry breaking was experimentally observed in the upper part of the 
transitional regime of pCf [28] , a range of Reynolds numbers which is out of reach of our 
present modeling due to its reduced cross-flow resolution. 




o = a2(/3-/3')f/iK-7oi?-'f/o. 



(26) 
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Wf) retain their initial values which, in turn, can freely be set to zero due to 
the in-plane Galilean invariance specific to stress-free boundary conditions. 

3 Numerical implementation 

The equations derived above have been implemented on domains of various 
streamwise and spanwise sizes and with periodic boundary conditions, 
which allowed us to develop a standard Fourier pseudo-spectral numerical 
scheme |29]. The diagonal viscous terms were integrated in spectral space 
exactly. The nonlinear advective terms and the non-normal terms arising 
from linearization around the base fiow were evaluated in physical space and 
integrated in time using a second order Adams-Bashforth scheme also used 
to treat the average velocity components Uj,Wj, j = 0,1. Fast Fourier 
transforms were used to pass from spectral to physical space and vice versa. 
Convergence of the results was checked by taking different time steps 6t and 
space steps 6x, 6z. In the range of Reynolds numbers under consideration, 
6t = 0.01 and 6x = 6z = 0.25 were found appropriate and used throughout 
the study. 

Simulations were performed in boxes with sizes ranging from V = x 

= 8 X 8 up to 64 X 64, although the model was originally designed for 
wide domains. A few short control runs were performed for V up to 128 x 
128. At any rate, results presented here are mostly dedicated to checking 
that the no-slip model captures the essential features of the globally sub- 
critical nature of the transition to turbulence in pCf, as a prerequisite to 
future work. In fact, the sizes considered here are already several times 
larger than the expected streak spacing ~ 6 [19]. In the terms of chaos 
theory, we would speak of "weak confinement" (see [2Sl Chap. 8] and below 
§4.21) in contrast with "strong confinement" legitimating the reduction to low 
dimensional differential systems. 

Obtained with limited computing power, simulation results presented here 
have already their own interest as will be shown in §4.41 but let us emphasize 
the fact that accepting low cross-stream resolution was the price to be paid for 
enlarging the domains up to sizes that are much larger than what can be dealt 
with in direct simulations with highly resolved cross-stream dependenceH 

Two further points are worth mentioning about the implementation of the 

■^In their study of the relaxation of obUque turbulent bands observed in the experiments 
[25] , Barkley and Tuckerman [3U] performed direct simulations in a quasi-one-dimensional 
box making an angle with the streamwise direction, very long in the direction perpendicular 
to the bands but a few MFUs long in the direction parallel to them, which represents an 
intermediate stage. 
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model. First, the continuity conditions fll3lll4p are dealt with by introducing 
appropriate stream functions ^'o; and velocity potentials $i such that: 

f/o = f/o-5.^o, 1^0=^0 + 9,^0, (27) 

and 

f/i = t/i + - 9,^1 , Wi=Wi + + , (28) 

with 

l3Vi = A<l>i . (29) 

Though not necessary on general grounds, the introduction of the additional 
quantities Uq, Wq, Ui, and Wi, is forced by our choice of in-plane periodic 
boundary conditions for the fields ^'q, ^'i, and $i. This is because these 
uniform components are generated by contributions to the potential $i and 
stream-functions ^'o, that would vary linearly with x and z, a spatial 
dependence that is precluded by their representations as Fourier series ex- 
pansions inherent in our computational approach. 

The velocity contributions derived from these equations are next reintro- 
duced in the expression of the nonlinearities when needed. The equations 
governing \l/o and are obtained by cross-differentiating and subtracting 
the equations for the velocity components in the usual way. Taking the di- 
vergence of these equations yields an equation for the pressure which is next 
used to determine the potential part of the velocity field accounted for by 
$1. The main advantage of introducing the fields \E'o, ^'i, and $i is that one 
can construct arbitrary but physically relevant divergence free velocity fields 
as initial conditions since this characteristic is built in their definitions. The 
uniform corrections {U i, Uq, Wi, Wq) are computed in parallel by integrating 
fl25|26l) and their W counterparts. The equations that have been numerically 
implemented are gathered in the appendix. 

4 Results 

4.1 Global sub-criticality 

Global sub-criticality of the model is depicted in Figure [2] which displays the 
behavior of the mean perturbation energy per unit surface as a function of 
the Reynolds number R. A discontinuous transition can indeed be observed 
around Rg ~ 175. Whether Rg represents a real threshold (called global 
stability threshold in previous studies, hence the subscript 'g') or just marks 
some crossover is presently a debated topic for plane Couette flow [HH 132] 
as well as for Poiseuille flow in a circular pipe which is also a linearly stable 
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Figure 2: Variation of the mean perturbation energy per unit surface with 
R. Evidence of a global stability threshold Rg ~ 175. Values corresponding 
to squares R > Rg correspond to sustained regimes, those corresponding 
to transients are marked with asterisks. Up and down triangles indicate 
the root-mean-square amplitudes of instantaneous fluctuations around the 
corresponding means. The computational domain is "D = 32 x 32. 



flow that becomes turbulent under finite amplitude perturbations [32l |33l 
IM] . For R > Rg spontaneous collapse of the uniformly turbulent state was 
not observed, even when pursuing the simulation for very long durations. 
For example, at i? = 175 the turbulent state persisted over more than 3 x 
10^ time unitsfl On the other hand, below ~ 174.5, the laminar state was 
systematically obtained but at the end of turbulent transients that could be of 
variable lengths. This result was obtained by combining quench experiments 
in which a state prepared a.t R = Ri = 200 was allowed to evolve after R had 
been suddenly decreased to some lower value i?f and adiabatic experiments 
where R was varied by small increments and the initial state at the new lower 
value of R was taken as the final state of the experiment at the previous value. 
The persistent turbulent regime at R = 175 was obtained in this second way, 
which — according to us — indicates a lower stability bound for the turbulent 
state at least over the corresponding time interval and for the considered size 
(Lj- = Lz = 32) since the procedure minimizes the perturbation brought to 

^For experiments in water with a gap 2h ~ Tmm, in the transitional range, R ^ 315, 
the time unit (turnover time h/Up) is 3.8 x 10~^s, so that the time length ^ 3 x 10^ would 
correspond to about 3 hours, much longer than what could be done in the laboratory 
under sufficiently well controlled conditions [55] . 
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the flow, especially regarding its uniform component Ui. In contrast, above 
Rg the initial finite perturbation brought to the system in quench experiments 
may be sufficient to make it leave the attraction basin of the turbulent state 
and return to the laminar state. Indeed, if the turbulent regime is a genuine 
attractor, its attraction basin shrinks to zero as R approaches -Rg from above 
and, with it, the size of the perturbation brought to the turbulent state 
necessary to make it decay toward the laminar state, in particular that linked 
to the fact that Ui is not at equilibrium for the new value of R after the 
quench. At any rate, it appears that the range of interest for the transition 
is no longer as low as in the sf case but approaches the experimental value 

A series of transients obtained in quench experiments from R[ = 200 
to -Rf = 171 is illustrated in Figure [31 It can be seen that they end quite 
abruptly so that it makes sense to define a conditional mean perturbation 
energy restricted to the plateau value, before final decayEI Such conditional 
averages are indicated as asterisks in Figure [2] while squares denote averages 
corresponding to sustained turbulent regimes. The line through the data 
points is a fit against a parabolic expression that has no theoretical founda- 
tion but strongly suggests that it goes smoothly through the critical value 
Rg. In particular, taking the mean perturbation energy as an order param- 
eter, we find no evidence of the proximity of a turning-point bifurcation as 
was tentatively suggested in earlier studies from analogies with elementary 
bifurcation theory [35J. 

The amplitudes of the fiuctuations of the average perturbation energy 
(conditional for transient states) are indicated by the up and down triangles 
in Figure [2] and, again, no remarkable behavior (i.e. divergence) can be 
noticed around i? ~ i?g. As discussed in the next sub-section, their amplitude 
is comparatively large only because, in this study, the size of the domain V 
has been chosen relatively small. The detailed statistics of the transients' 
lifetimes, which is the main original result presented in this paper, is studied 
in gai 

4.2 Extensivity of the sustained turbulent regime 

It has been argued elsewhere by one of us (PM in [2j) that the modeling in 
terms of low- dimensional ordinary differential systems [12] niay give valuable 
hints only about the mechanisms of turbulence sustenance well beyond the 
transitional regime but not necessarily about the transitional regime itself, 

similar observation was made by Faisst & Eckhardt '.'At, in numerical simulations 
of Poiseuille pipe flow. They pointed out that it was not possible to distinguish between 
sustained turbulence and turbulent transients before decay. 
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Figure 3: Total perturbation energy per unit surface Etot as a function of 
time in a series of transients from a turbulent flow prepared at R[ = 200 
suddenly quenched at -Rf = 171 < -Rg ~ 175. 



Table 2: Perturbation energy per unit surface: time average and fluctuations 
for R = 200. 
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where spatio-temporal behavior is involved. Indeed, the main underlying 
hypothesis of this kind of modeling is that the space-time dependence of the 
perturbations can be described by mixing a small number of modes with 
frozen space dependence and time varying amplitudes. Specific resonances 
between modes make the temporal properties of the system, notably the 
fractal properties associated to chaotic transients, excessively sensitive to the 
physical size of the equivalent confined system. Furthermore, the approach is 
intrinsically unable to deal with the growth or decay of turbulence through 
the coexistence of laminar and turbulent domains that fluctuate in space 
and time [40J. When considering this speciflc problem, it seems legitimate 
to require that, to be appropriate, a model should display some kind of 
statistical robustness when the size of the simulation domain is varied [3H] . 
An indication that our model behaves appropriately in the turbulent regime 
is obtained by considering the total perturbation energy per unit surface and 
its fluctuations as the surface of the simulation domain is increased. Table [2] 
displays our results. A series of experiments over domains V = L^. x 
was performed, with domain surfaces multiphed by factors n = 1, 2,. . . , 64, 
starting from P = 8 x 8 and for R = 200. 

For the smallest systems considered, turbulence was only transient at 
R = 200. Needing larger R to get sustained turbulence is not surprising since 
much dissipation is associated with the in-plane space dependence forced by 
periodic boundary conditions at short distances, which raises the thresholds 
for complex behavior. For P = 8 x 8, the domain can only flt a very small 
number of structures so that the transient is better analysed in terms of 
chaotic saddles in the context of temporal chaos [39] . For D = 16 x 8 turbu- 
lence is again transient but the dynamics has a more spatio-temporal flavor 
than for 8x8, hence closer to what is obtained at larger sizes. Evidence that 
the turbulent regime at i? = 200 is extensive for larger domains, in practice 
beyond P = 16 x 16, is obtained from the facts that (^) the time average of 
Etot is independent of the size and (ii) the fluctuations around this average 
value decrease roughly as the inverse square root of the size. 

The extensivity property is valid for the uniformly turbulent regime but 
size effects may still be expected to affect the transition process, especially in 
view of the presence of large scale flows resulting from the space modulation 
of the mean flow correction to be discussed in the next sub-section, as soon 
as strong inhomogeneities are present. This study only suggests us a typical 
size below which the temporal approach is better suited. 
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Figure 4: Correction to the base flow for a sustained turbulent regime ob- 
tained as described in the text {R = 200, V = 32x 32). Top panel: Traces of 
the average flow components Ui, Wi, Uq, and Wq as functions of time. Bot- 
tom panel: Mean turbulent flow profile = y + UiCy{l — y"^), the value of 
Ui is obtained by averaging the above relevant trace over t G [0.5, 15] x 10^. 

4.3 Mean flow in the sustained turbulent regime 

An important feature of our model is that it already contains corrections to 
the base flow at the lowest truncation order. Figure H] illustrates the result of 
a specially designed experiment at i? = 200 starting from an initial condition 
taken from some fully turbulent state, in which all the components of mode '1' 
were artificially turned off. During a brief transient that has been cut off, the 
streamwise contribution to the mean flow Ui builds up. As seen in the figure, 
it is negative and statistically constant, with time average {Ui)t = —0.221 
and standard deviation Ejj^ = 0.015. In contrast, the transverse contribution 
averages to zero: {Wi)t = 1.4 x 10"^ with = 6.9 x 10"^ 50 {Wi)^ 
for the specific numerical simulation considered). At the same time, Uq and 
Wq also average to zero: {Uo)t = -2.4 x 10"^, S^^ = 3.3 x 10"^ and Wq = 
—1.5 X 10~^, ~ ^-^ ^ 10"'^. The mean turbulent flow profile is displayed 
in Figure H] (bottom) as the superposition of a correction Cy{l — y) with 
amplitude Ui to the base flow Uh = y, pointing out the expected formation 
of a central region with reduced shear. It will be interesting to study the 
patches of negative Ui that appear during the growth/decay of turbulent 
spots [37], as a local counter-part to the steady state uniform correction in 
the sustained turbulent regime. 
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4.4 Transients below the global stability threshold 

The global stability threshold -Rg is best defined as the value of the control 
parameter below which the nontrivial state, here the turbulent regime, is 
unstable and thus cannot be sustained in the long term. Accordingly the 
trivial state, here the laminar fiow, is the only possible equilibrium solution, 
whatever the initial configuration (or accidental perturbations). Many stud- 
ies have been devoted to the transition from laminar to turbulent fiow, either 
under the effect of natural fluctuations, upon triggering speciflc localized 
perturbations, or upon modifying the base flow in a well-controlled way (see 
the review in plj). In such experiments, the turbulent state may be reached 
with flnite probability only due to the local stability of the base flow, which 
explains the dispersion of results in early studies and some sensitivity to the 
triggering process. In contrast, from the very deflnition of Rg, it may seem 
less ambiguous to start from the turbulent side at large Reynolds numbers 
and study its relaxation as R is decreased. As already mentioned one can 
consider either an adiabatic decrease of R which minimizes the risk of intro- 
ducing dangerous perturbations, or else quench experiments using the same 
methodology as the one which produced the experimental value Rg ~ 325 
|19j . The cumulative distribution of transient lifetimes was then studied and 
found to display an exponential tail of the form n(r' > r) oc exp(— r/r/j), 
where n(r' > r) is the probability of observing a transient with length r' 
larger than some given time r. In this expression tr is the characteristic de- 
cay time of the distribution. If it was effectively exponential from the start, 
then tr would be the mean transient length at the corresponding value of 
R. From that estimate, it was proposed that tr diverges to infinity as R 
approached -Rg from below as 

(n,)~(i?g-i?)-i, (30) 

see [35l [3T] . This behavior was recently questioned and a reexamination of 
the data lead to the proposal of an indefinite exponential increase [32j though 
the time scale computed in notelUleaves little hope to check the extrapolation 
at larger R experimentally. Our quench experiments with the ns model is a 
tentative contribution to the debate. 

Experiments were performed by quenching the system from series of uni- 
formly turbulent states at R = Ri = 200. These states were obtained from 
snapshots taken during a single long experiment and sufficiently separated 
in time to be considered as independent turbulent initial conditions. The 
Reynolds number was next suddenly decreased to a final value R = Rf. In 
order to be able to perform statistics over a large number of independent tri- 
als, we took a domain size T> that was reduced as much as possible, though 
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Figure 5: Distribution of transient lifetimes in quench experiments with = 
200 and variable Rf in semi- log plot. 

sufficiently large to preserve the extensivity property of the uniformly turbu- 
lent regime. P = 32 x 32 was chosen but runs with V = 128 x 32 or 128 x 64 
were performed for purposes of control or pattern imaging (see §4.51) . 

The lifetimes of individual transients were determined as the times needed 
to reach some cut-off value of the total perturbation energy per unit surface 
Et below which the system immediately and irreversibly decays to the lami- 
nar state, as can easily seen in Figure [H] displaying a typical set of experiments 
to be analysed. The statistics were checked to be insensitive to the precise 
value chosen for the cut-off provided that it was sufficiently below 0.015. 
Here the cut-off was uniformly taken equal to 0.0015 (see also later, Fig. [TOl) . 
The very late stage of decay for Etot < 0.0015 appears to correspond to a 
fast exponential viscous damping of streaks. 

Below R = 175, transients with variable lifetimes began to appear and 
well-defined lifetime distributions were seen to develop, while relaxation could 
be considered as "immediate" below i? = -Ru ~ 164 (notation introduced in 
|35j). Upon approaching Rg ^ 175 from below, using lin-log scales Figure [5] 
shows that the statistics of transient lifetimes are indeed exponentially de- 
creasing and that the slope decreases as Rf increases. 

The slopes of the exponential tails were determined by fitting them with 
straight lines. The result is plotted in Figure [6] as open dots, either in linear 
scale (top) or semi- log scale (bottom). If the distribution were exponential 
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Figure 6: Slope of histogram tails in lin-log scale (Fig. E]) as a function of R. 
Top: lin-lin. Bottom: lin- (natural) log scale. Open dots: slope of histogram 
tails in semi- log plot; open squares: inverse mean length of transients com- 
puted as indicated in the text. 

from the start, up to some offset value r^in, the slope would be given by 
l/(r — Tmin) where (. . .) is the arithmetic mean over the observed transient 
lengths r. Corresponding results are plotted as open squares in Fig. El Both 
estimates compare well with each other down to i? ~ 169 below which they 
diverge. In fact, histograms for R < 170 were not presented in Fig. in order 
not to overcrowd it with curves that would have needed another scale to 
make them discernable. These histograms display a shoulder corresponding 
to an accumulation of short transients that biases the computation of the 
slope from the inverse mean and explains the discrepancy between the two 
estimates. Conversely this assures us that the fit to an exponential is indeed 
satisfactory for R > 170. 

It is, however, hard to put error bars on our results because we have a fi- 
nite number of transients at our disposal, usually between 100 and 200 except 
for R = 174 and 174.5 for which we have only about 30 transients. To un- 
derstand the origin of the difficulties, one must consider how the histograms 
evolve as the number of experiments increases. For example, the presence 
of a small number of measurements that seem atypical of the statistical tail 
implies an increase of the number of experiments by a large factor to rub out 
their effect. All the values plotted in Figure [6] are what we consider to be 
our best estimates based on careful inspection of individual histograms and 
comparisons of statistics obtained by including or excluding transient lengths 
that might be judged atypical of the overall decay behavior. In general all re- 
sults were consistent to within less than 10%, especially at the largest values 
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Figure 7: Slope histogram tails in lin-log scale as a function of R. Enlarge- 
ment of Fig. [HI Top: lin-lin scale. Bottom: lin-(natural)log scale; the added 
dashed line is a guide for the eyes, showing the deviation of the top highest 
points from the overall exponential behavior. Symbols as in Fig. El 

of R considered. 

This being said, the zoom on the variation of the slope as a function of R 
displayed in Figure [7] confirms the exponential decrease of the terminal slope 
of histograms, or conversely an exponential increase of the characteristic de- 
cay time. However, the values corresponding to the points at i? = 174 and 
174.5 seem somewhat too small for them to stay aligned with the other points, 
as suggested by the dashed line in Figure [7] (bottom). This line is obtained 
from a linear fit of the logarithms of the slopes measured for i? G [170, 173.5], 
which we consider as our most reliable estimates. The deviation can have 
two explanations: 

- either, as discussed earlier, it is an artifact of our estimation procedure 
with an insufficient number of points but it turns out that the corresponding 
distributions were surprisingly regular and the slope estimation unambigu- 
ous, 

- or it is a real effect indicating that the characteristic time increases faster 
than expected, i.e. diverges for some finite -Rg that should be close to 175 
(for which the turbulent attractor was found to be stable). 

For R = 174 and 174.5, an over-representation of long transients among 
our 30 or so samples would indeed explain the observed under-evaluation 
of the slope, which might be the trace of a genuine (turbulent) attractor for 
slightly larger values of R. Note also that the difference between the estimates 
obtained by slope fitting and by inverse means at i? = 174 — quite anomalous 
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Figure 8: Mean perturbation energy as a function of time during a transient 
in a quench experiment at R = R{ = 168 starting from Ri = 200 for = 128, 
= 64. 

in view of the general agreement between the two estimates when R > 170 — 
is entirely accounted for by the accidental absence of a sufficiently short 
transient Tj^m- Accordingly, we do not exclude the possibility of a crossover 
between the exponential regime typical of transients tentatively attributed 
to a chaotic saddle and a real divergence at some well defined Rg above 
which a turbulent attractor does exist. This divergence would be typical of 
a crisis converting the chaotic repellor into an attractor. The anomalously 
long transients observed for R = 174 and 174.5 would then be representative 
of the behavior just below the crisis point. 

Finite size effects, still sizable at "D = 32 x 32, will make it difficult 
to decide what is the actual behavior but we believe that considering larger 
systems, which is presently in progress, will help us to obtain a clearer picture 
of the transition. The next section is thus devoted to an illustration of 
a transient in a system of intermediate size with the purpose of showing 
that one cannot avoid the spatio-temporal feature of the dynamics in the 
transitional range of Reynolds numbers. 

4.5 Illustration of decay during a transient at R = 168 

In this subsection we consider a typical transient obtained for R = 168 with 
Lj. = 128, = 64, starting from R = Ri = 200. Figure [H] displays the 
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Figure 9: Four snapshots taken during a transient in a domain V = 128 x 64. 
Dark/deep red (clear/light yellow) regions correspond to small (large) pertur- 
bation speed [colored version on line]. One point out of 4 in each direction 
is represented, hence Ax = Az = 1. The horizontal axis corresponds to 
the streamwise direction. The contribution of streaks to the perturbation is 
clearly visible at the laminar /turbulent frontier. 

evolution of the mean perturbation energy as a function of time throughout 
the transient. The general aspect of this curve is similar to those in Figure [3], 
with a distinct plateau followed by an irreversible decay. Here, the transient 
is not very long since R{ is not close enough to Rg. Snapshots of the local 
perturbation velocity amplitude u{x, z, t) are given in figure [9l This quantity 
is defined as u{x, z, t) = ^2{Eo{x,z,t) + Ei{x,z,t)), where Eq = 1{U^ + W^) 
is the local kinetic energy contained in the streaks and Ei = ^{Ui +Vi +Wi) 
the kinetic energy contained in the remaining part of the departure from 
the laminar profile, including the contribution of streamwise vortices. This 
definition is reminiscent of that of a turbulent velocity intensity except that 
the latter is defined with respect to the mean flow. Here, this would have 
led us to subtract Ui, but would no longer have made sense for the spatially 
inhomogeneous turbulent regime during decay, as illustrated below. 

Figure M presents several snapshots of the perturbation flow pattern at 
different key instants. First for t = 100, just after the quench, turbulence 
initially prepared at = 200 gets stabilized at the level corresponding to 
-Rf = 168. Though it is relatively uniform, one can already detect some 
large scale inhomogeneity. Later, this inhomogeneity increases up to a point 
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Figure 10: Variation of the total perturbation energy during the end of 15 of 
the transients in Fig. [SI 

when obhque bands are conspicuous, as in the second image at t = 800 
while the level of perturbation energy has remained essentially the same. 
Subsequently, the turbulent region becomes narrower and the perturbation 
energy decreases accordingly. At t = 1500, is about the half of the initial 
energy content. The regular regression of the band ends when it breaks into 
several parts forming smaller patches that, in turn, decay as seen here for 
t = 2000. The regular regression of the width of the turbulent bands explains 
the roughly linear decrease of the mean perturbation energy for T > 1000 in 
Fig. [H] assuming that the front between the turbulent and laminar domains 
moves at constant speed. 

In fact, a similar behavior was probably exhibited by transients already 
for P = 32 X 32 though the picture in physical space is less clear. Indirect 
evidence can be obtained by zooming on the last 500 times units of each 
transient. Figure [TU] displays a superposition of 15 of them for R = 171 
taken from the series used for Fig. [SI Statistically speaking, they all follow 
the same roughly linear decay, from the end of the plateau to the cut-off cor- 
responding to the beginning of the final exponential viscous damping (not 
represented). Accordingly, and comparing to what happens for V = 128 x 64, 
it seems that we have a nucleation problem with a probability distribution 
for laminar germs of different sizes so that when a wide enough germ appears, 
the energy plateau ends and the turbulent domain recedes regularly. (At this 
stage, it could be noticed that uniformly subtracting about 500 time units 
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from the transient lengths would not change the overall statistical behavior.) 
The intensity of the mechanism sustaining turbulence, a function of R, is 
certainly essential in controlling this probability. Such a mechanism is some- 
times studied using low dimensional dynamical systems modeling p^fTH fTG]. 
Accordingly, a way to reconcile this picture with that of chaotic saddles could 
be recovered but it is still not clear how to reintroduce physical space and 
correlative size effects. At any rate, considering a larger system will certainly 
shed light on this nucleation probability problem. 

It is worth to mention that, during the decay and even at rather late 
stages, the energy remains much concentrated in the turbulent patches. It 
was empirically found that setting the cut-off between turbulent and laminar 
flow at M = O.SMref, whcrc Mref is the mean of u in the starting reference 
state at t = 100, well discriminates the two regimes from each other. For 
example, at t = 2000, one finds that about 70% of the remaining energy is 
concentrated in 10% of the total surface. Going backwards in time, one gets 
that, at t = 1500, the turbulent band occupies about 48% of the surface and 
contains around 90% of the energy, and still earlier, that for t = 800, the 
turbulent fraction is about 87% and contains 98% of the energy. 

To conclude this section, it should be remarked that the oblique band 
observed in Fig. [9] has no connection with the oblique turbulent bands ex- 
perimentally observed in the upper transitional range ^28j and numerically 
reproduced by Barkley and Tuckerman [SO]- Here, its presence is clearly 
linked to the periodic boundary used and its development directly follows 
from the oblique ovoid shape of the inhomogeneity already present at the 
very beginning of the quench experiment. Depending on the size and aver- 
age orientation of this inhomogeneity, the ultimate shape of the turbulent 
domain is either a longitudinal, a transverse, or an oblique band. In view of 
the size of modulations at t = 100 (Fig. [HI top), we think that a system at 
least four times as large as that used for this illustration should be considered 
in order to have a chance to explore several inhomogeneity patterns before 
decay. This is what motivates the simulations over a domain T> = 256 x 128 
presently in progress and alluded to above. 

5 Conclusion 

In this paper we have considered the low-i? range of transitional pCf. On the 
one hand, modeling this situation in terms of low dimensional deterministic 
dynamical systems is limited due to the inherently spatiotemporal charac- 
ter of the problem. For example, decay from the sustained turbulent state 
close to the global stability threshold Rg is not well described as a chaotic 
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transient close to a crisis bifurcation point, even though exponential distri- 
butions of transient lifetimes are observed near this point [T3], because this 
concept introduced in nonlinear dynamics is not adapted to the account of 
decay through regression of turbulent patches fluctuating in time and space. 
On the other hand, the direct connection to statistical physics suggested by 
Pomeau ^40j and underlying the abstract spatio-temporal intermittency ap- 
proach [H] may seem far-off though it takes spatial extension into account 
in an analogous way. Extending previous work by one of us [20l [22], the 
approach developed here intends to bridge the gap between low and high di- 
mensional systems in a concrete way and with a semi- quantitative ambition. 

Except for the fact that it fell short of predicting reasonable values of the 
transitional Reynolds numbers, the early stress-free attempt already con- 
tained many interesting features and its deficiency could easily be attributed 
to the use of unrealistic boundary conditions in the derivation of the model. 
The first part of the present paper was thus devoted to the derivation of a 
similar model, but appropriate to the physically realistic no-slip boundaries 
[211 [2Z! • When truncated at the lowest significant order ( §2.21) . this new 
model was seen to have nonlinear and non-normal structures close to those of 
the stress-free model but with slightly different coefficients when basis func- 
tions used in the Galerkin expansion were appropriately normalized. The 
major differences, discussed in §2.31 appeared to be (i) the presence of an 
additional damping of the drift flow that was no longer cross-stream indepen- 
dent but had a more generic parabolic profile, (m) a strongly enhanced viscous 
dissipation of the other flow components (associated to steeper variations of 
the Galerkin basis functions close to the plates), and {Hi) the existence of a 
non-trivial mean correction to the base flow Ui already at lowest truncation 
order. As could be expected from the structural similarity of the models, the 
same qualitative behavior was obtained but a wide step in the right direction 
was made since the study presented in the second part of the paper showed 
that transitional Reynolds numbers of interest were now less than a factor of 
two off the experimental range, instead of a factor for ten with the stress-free 
model. 

Experiments reported in [121112] clearly showed that, in the low-i? part of 
the bifurcation diagram, the dynamics is dominated by processes involving 
flow structures that are coherent over the full gap between the plates and 
this is certainly the reason why the results obtained from the model are 
reasonably satisfactory. Increased dissipation associated to more resolved 
cross-stream dependence is expected to improve the agreement at the price 
of much more analytical and computational work already around Rg. 

On the other hand, increasing the order of the model would certainly be 
necessary much beyond Rg, in particular in the range where the transition 
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'oblique turbulent bands ^ featureless turbulence' takes place [SO]- This 
is however not advisable since the model would then get more and more 
cumbersome explicit nonlinear terms from the Galerkin expansion which, 
in turn, would be truncated at increasingly insufficient orders, while it is 
clearly more economical to use a high precision pseudo-spectral scheme in 
the cross-stream direction. 

Besides the observation that fully developed spatiotemporal chaos is ex- 
tensive ( §4.21) . the main concrete result presented here relates to the determi- 
nation of a possible global stability threshold Rg and the statistical behavior 
of transients close to Rg for systems of moderate size. As in laboratory 
experiments, a discontinuous transition was observed ( §4.11 Fig. [2]). Tur- 
bulent transients accounting for the relaxation of turbulence during quench 
experiments were shown to exhibit exponentially distributed lifetimes ( §4.41 
Fig. [5]). Beyond R ~ 175 the turbulent state was observed over time inter- 
vals that could seemingly be arbitrarily long, hence suggesting Rg = 175. For 
R < 175, the slopes of the exponential tails of the lifetime distributions were 
seen to approach zero exponentially as R increased, with a possible crossover 
to a diverging behavior close to the supposed Rg, in spite of post-processing 
uncertainties (Fig. [7j). 

Extending this conclusion to the case of transitional pipe Poiseuille flow 
[321 [331 IM] needs additional investigation before tentatively establishing any 
connection. Such a connection, which should be taken with care, has been 
made from a reanalysis of the pCf results that also leads to indefinite expo- 
nential decay with R of the slopes of histograms' tails. In our opinion, in the 
Couette case available experimental data were however obtained for values of 
R somewhat too far from the hypothesized Rg ~ 325 (determined from sev- 
eral independent ways) for supporting the conjecture of an underlying chaotic 
saddle typical of temporal chaos, especially in view of experimental pictures 
[35] . On the contrary, in close correspondence with these pictures, the il- 
lustration given in §4.51 clearly points to spatio-temporal chaos as explaining 
the transient behavior, with the statistical possibility of a local breakdown 
of the turbulent state, further extending to the whole system. This behav- 
ior is reminiscent of what happens at a thermodynamic phase change [10] . 
The size of the system for which these preliminary results have been ob- 
tained, V = 32 X 32, though situating it apparently inside the extensive 
regime, is clearly insufficient to probe any thermodynamic-like property of 
the transition, and a fortiori any issue related to universality in the sense 
of critical phenomena in statistical physics. We expect that the study cur- 
rently in progress for V = 256 x 128, a size unreachable by direct numerical 
simulations, will bring us interesting results. 

Since it represents a simplified version of the Navier-Stokes equations 
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at least for R not too large, our model can serve to study other questions 
posed by the transition to turbulence in globally sub-critical flows. For ex- 
ample, whereas it is clear how streamwise vortices generate streaks through 
the lift up mechanism introduced long ago, instability mechanisms and some 
processes involved in the regeneration of vortices are still unclear in spite 
of recent progress P|. The very coexistence of laminar and turbulent flow 
implied by global sub-criticality, the mechanisms by which one of the phases 
gains against the other, i.e. how spots grow {laminar —>■ turbulent transi- 
tion) or how the turbulent state recedes {turbulent — >• laminar transition), 
how asymmetric these problems are, what is the role of the base flow modi- 
fication inside a turbulent domain (local counter-part to §4.31 in the presence 
of large scale modulations), are questions that we intend to examine, with 
the stimulating perspective of going beyond the case of pCf and considering 
other wall flows of less academic interest but experiencing a similar transition 
to turbulence. 



Appendix: Equations used in simulations 



The equations for \E'o, \E'i, and $i that have been effectively implemented 
read: 



dt-R~'{A-jo)] ^0 

dt-R-'{^-li)] ^1 
dt - R-^ (a - p^y^ (a^ - 2/3^ A + 7i/3' 



$1 



-1 r 



A -13'] A-'l3\d,Nu, + d^Nw,) - l3Ny, 



where (3 = = y/S everywhere and the expressions of the advection terms 
have been rewritten in the energy conserving form mentioned in the text: 



Nwo 
Nw, 



Mv, + (3E^ , 
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where the "weighted" energies read: 

= \ai{Ul + W^) + \a2{Ul + Wl) + \a:iV^, 
= a2(Wi + H^oW^i) , 

and 

Mwo = aMd,Wo-dM+c^2U,id^W,-d,Ui) 

+ Vi{a2P'W,-a,d,Vi), 
Mu, = a2[Wo{d,U,-d,W,) + W,idM-d^Wo)] 

- a2P"ViUo , 

Mw, = a2[Uoid,W,-dM,) + U,idM-dM] 

Mv, = U^{a;d^Vi-a2l3Ui) 

+ Wo{a^d,Vr - a2(3Wi) ■ 

Velocity components are given by (!271 - [29i) . 

The evolution is computed in Fourier space by a numerical scheme that 
treats linear terms exactly and nonlinear terms by an Adams-Bashforth for- 
mula. Formally writing the evolution problem as 

f^X = CX + X{X) 

it can be seen that, in order to keep second order accuracy, we have to 
compute 

X{t + 5t) = exp{C5t) X 

[X{t) + 6t{1.5Af{X{t)) - 0.5exp{LSt)Af{X{t - 6t)))]. 

When integrating the equations, it can also be verified that their r.h.s. cancel 
exactly at wave-vector (fc^, k^) = (0, 0) so that taking the inverse of A is a 
legitimate operation. As stated in the main text, average velocity compo- 
nents must be treated separately. Corresponding formulas for {Ui,Uo) are 
obtained from ( 125)1261) and those for and {Wi, Wq) by two similar equations 
with U replaced by W. 
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